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Abstract. An accurate MonteCarlo simulation of 
the deficit of primary cosmic rays in the direction 
of the Moon has been developed to interpret the 
observations reported in the TeV energy region until 
now. Primary particles are propagated trough the 
geomagnetic field in the Earth-Moon system. The 
algorithm is described and the contributions of the 
detector resolution and of the geomagnetic field are 
disentangled. 
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I. Introduction 

Since the galactic cosmic rays are hampered by the 
Moon, a deficit of cosmic rays in its direction is expected 
(the so-called "Moon shadow"). The Moon shadow 
is an important tool to calibrate the performance of 
an air shower array. In fact, the size of the deficit 
allows a measurement of the angular resolution and 
the position of the deficit allows the evaluation of the 
absolute pointing accuracy of the detector. In addition, 
charged particles are deflected by the geomagnetic field 
by an angle inversely proportional to their energy. The 
observation of such a displacement provides a direct 
check of the relation between shower size and primary 
energy, thus calibrating the detector. 

A detailed Monte Carlo simulation of cosmic ray 
propagation in the Earth-Moon system is mandatory to 
understand the Moon shadow phenomenology and to 
compare the observed westward displacement with the 
expectations in order to disentangle the geomagnetic 
effect from some possible experimental biases. 

II. Monte Carlo Simulation 

A. Simulation strategy 

When the emission of photons by a given gamma- 
ray source has to be simulated, a classical approach can 
be followed: once its energy spectrum is known, it is 
enough to decide for how long the source emission has to 
be reproduced to calculate the number of expected events 
which have to be sampled. Of course, if the response of 
an Earth-based array has to be simulated, for each one of 
these primary particles also the induced shower and the 
detector response have to be calculated. The easiest way 
to do that, is to shoot up the particles from the detector 
to the source. This trick allows to concentrate the efforts 

'Actually, from a suitably-choosen area around it. 



of the simulation only on the showers which effectively 
construct the signal. 

Nonetheless, this method cannot be used for the 
Moon-shadow as it is, because of two main reasons. 

1) The Moon shadow is not a signal, but a "lack" in 
the background. 

2) the effect is provoked by charged particles and not 
by photons. This implies that we must take into 
account the effect of the electro-magnetic fields 
on their trajectories. 

The first argument is the more relevant one. In fact, 
if we wanted to remain faithful to the approach of 
simulating only the particles which reach the detector, 
we should compute the showers of the background 
sorrounding the Moon and then take off the ones coming 
exactly from within Moon. In such terms, the simulation 
is likely unfeasible, at least for low energy threshold ex- 
periments. To be precise,if the detector energy threshold 
is a few hundreds of GeV, even by considering a square 
sky window sorrounding the Moon not so large (e.g. 
10°), too many showers should have to be simulated to 
reproduce even just one year of data taking. 

It is evident that the Moon shadow requires a different 
strategy of simulation. It is better to treat the Moon 
like a standard source and then reverse the amplitude 
of the signal. In the end, both the gamma-source case 
and the Moon shadow case reduce to a perturbation of 
the cosmic -rays background. The only difference is the 
sign of the perturbation. 

Because of their electric charge the cosmic rays do 
not proceed straight from the Moon to the Earth, but 
their trajectories are bent by the geomagnetic field. Since 
this effect plays a crucial role in the final result of the 
physics analysis, it must not be underestimated during 
the simulation. Especially for those particle which have 
low energy, the bending effect is very strong and a 
realistic prediction is possible only if the simulation 
accounts correctly for the intensity and the direction of 
the geomagnetic field. What is more, notice that when 
the particles are sent back to the Moon from the detector, 
the sign of their charge must be reversed to properly 
reproduce the direction of the deviation. 

B. The generation and the detection of the showers 

The air showers development in the atmosphere has 
been generated with the CORSIKA v. 6.500 code in- 
cluding the QGSJET-II.03 hadronic interaction model 
for primary energy above 80 GeV and the FLUKA 
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code for lower energies [1]. Cosmic ray spectra have 
been simulated in the energy range from 30 GeV to 1 
PeV following the relative normalization given in fl2), 
resulting from a global fit of main experimental data. 
About 10 8 showers have been sampled in the zenith 
angle interval 0-60 degrees. The secondary particles have 
been propagated down to a cut-off energy of 1 MeV. 

At present, the only experiment able to observe the 
Moon shadow with high statistical significance and 
energy threshold well below 1 TeV is ARGO-YBJ 0. 
Therefore we reproduced an ideal detector placed in 
YangBaJing (4300 m a.s.l.) having geometrical features 
similar to ARGO-YBJ and a duty-cycle of 90%. We used 
a GEANT4-based code [4|. The trigger threshold has 
been set to 20 charged particles over the whole detector. 
The core positions have been randomly sampled in 
an energy-dependent area large up to 10 3 x 10 3 m 2 , 
centered on the detector. 

C. The geomagnetic model 

It has been already noticed that if a primary cosmic 
ray (energy E, charge Z) traversing the geomagnetic 
field is observed by a detector placed in YangbaJing, 
its trajectory shows a deviation along the East-West 
direction[5||j which in first approximation can be written 
as: 



The sign is set according to the usual way to represent 
the East- West projection of the Moon maps (see Fig. 0). 
The eq. [T] can be easily derived by assuming that the 
geomagnetic field is provoked by a pure dipole laying 
in the centre of the Earth (see Appendix). Nonetheless 
this approximation, which is derived for nearly vertical 
primaries, is not enough when the primaries energy is 
below few TeV. 

To perform a numerical estimate of the bending effect, 
it is necessary to adopt a model of the magnetic field 
in the Earth-Moon system. The roughest one is the so- 
called Virtual Dipole Model (VDM), whose name is 
self-explaining. A better choice is the Tsyganenko-IGRF 
model (T-IGRF hereafter) [6|, which accounts for both 
internal and external magnetospheric sources. We com- 
pared the effect on the particle trajectories of VDM and 
T-IGRF, in both cases finding non negligible differences 
with respect to the -1.58°Z/E[TeV] formula. Among 
the two models themselves, we observed discrepancies 
up to the ~ 15% level (4° 7°) for sub-TeV primary 
energies, mainly due to the field intensity near the Earth 
surface. Since the T-IGRF model accounts for more 
factors, we refer to it hereafter. In Fig. Q] you can 
appreciate the analytical trend (eq. [TJ together with the 
actual East- West displacement calculated applying the T- 
IGRF model for protons and He nuclei. The analytical 
approach clearly underestimates the East- West deviation. 

2 No deviation is expected along the North-South one. 
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Fig. 1: Deviation induced by the geomagnetic field on 
protons and He nuclei. Each point refers to a simulated 
shower. Both the arrival direction and the date of the 
propagation are randomly sampled, respectively from an 
isotropic distribution within the sky and from a uniform 
distribution over 2008. The analytical trends obtained 
from the equation Q] are also shown. 




Fig. 2: Residual displacement with respect to the analyt- 
ical expectation. The deviation is calculated by applying 
the the T-IGRF model (see text). The color scale repre- 
sents the number of showers laying on the single pixel. 

The Fig. [2] shows the differences of the T-IGRF- 
induced deviation with respect to the leading term 
1.58 Z/E[TeW]. The upper (lower) panel contains such 
a residual deviation along the East- West (North-South) 
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direction as a function of the primary energy. Altough 
there are no effects for energies E > 10 TeV, below 
few TeV the residual displacement can reach 1°. Notice 
that unlike the analytical approach would suggest, the 
North-South deviation of the single primary is non-null, 
being zero only on average. 

D. Moon shadow simulation 

By following the procedure described before, we can 
obtain the Moon shadow maps represented in Fig. [3] 
There can be appreciated the displacement induced by 
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Fig. 3: Folding different contributions to the Moon 
signal. Upper part of the figure: Moon as it would be ob- 
served by a perfect detector without geomagnetic field; 
only the magnetic field is switched on. Lower part: only 
the detector is switched on; both the magnetic field and 
the detector are switched on. The maps are drawn using 
equatorial coordinates (declination VS corrected right 
ascension). Only the showers triggering the detector are 
shown (N > 20). The color scale represents the intensity 
of the signal. 



point-like source and its own finite angular width (half 
a degree) contributes to the spread. The superposition 
of the two effects can be easily visualized in case of 
gaussian Point Spread Function (PSF): 



RMS = 0-4/1 



0.13 c 



where the root mean square of the signal RMS is related 
to the variance a 2 of the PSF. The contribution of the 
Moon size to the RMS is (not) dominant when a is low 
(high), i.e. at high (low) particle multiplicities. Just to 
be explicit, the difference between RMS and a is 20% 
if (7 = 0.2°, less than 5% if a > 0.4°, and only 1.7% if 
cr = 0.7°. 
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Fig. 4: Effect of the PSF along the east-west direction. 
The continous black line represents the Moon shadow 
deformed by the geomagnetic field as it would appear 
to an ideal detector. By considering also the effect of 
the detector PSF, the diplacement of the signal peak 
is enhanced, moreover the well-known smearing effect. 
The figures represent only protons. 



the geomagnetic field. The long tail of the left part of 
the up-right map is made by the lowest energy particles 
(below 1 TeV) which are more deviated. Concerning 
the bottom-left map, the detector by itself provides 
the smearing of the signal, leaving intact the circular 
symmetry, as expected. 

It is possible to study the effects of the finite angular 
resolution of the detector and of the geomagnetic field 
separately. 

As already noticed, if we consider the magnetic de- 
viation but not the detector, the circular symmetry of 
the signal is broken only along the East- West direction 
(see Fig. [3] top-right map). That make us confident 
that only the smearing due to the angular resolution 
affects the signal along the north-south direction, thus 
allowing its determination. Actually, we stress that what 
we determine by considering the spread of the signal 
along the North-South direction, cannot be properly 
named angular resolution, because the Moon is not a 



10 



0.2- 



EW displacement 



-0.2 



-0.4 - 



I 



number of detected particles 

Fig. 5: Expected displacement of the Moon shadow in 
the East-West direction as a function of multiplicity. 
The upper scale refers to the median energy of rigidity 
(TeV/Z) in each multiplicity bin (shown by the horizon- 
tal errors). 

Finally, it is worth to invert the switching order of 
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the detector effect and of the geomagnetic field. In 
Fig. [4] the effect of the angular resolution on the East- 
West projection of the Moon is shown. Because of the 
signal asymmetry induced by the magnetic field, such an 
effect provides not only the smearing, but also a further 
displacement of the signal peak. The West tail of the 
shifted signal, in fact, has a larger weight than the sharp 
East one and tends to pull the signal in its direction. 

In Fig.|5]the total eastward displacement is plot versus 
the number of detected particles. 

III. Conclusions 

We reproduced the Moon shadow effect as observed 
by an ideal detector located at high altitude. Much 
attention has been devoted to the simulation of the 
geomagnetic field bending the particle trajectories. We 
have been able to disentangle the contribution to the final 
signal of the geomagnetic field and of the detector PSF 
respectively. We quantified the contribution of the finite 
Moon disc to the angular resolution estimation. For a 
comparison with the data collected by the ARGO-YBJ 
experiment, see 0. 

Appendix 

Here is shown how to obtain the formula[T] Since only 
the magnetic field is supposed to act upon the particles 
trajectories, the Lorentz equation can be written as: 



x(i) = x + v 1 



zee 2 r r dx 

I At j da — x B(x, a) 

E Jo Jo da 

(2) 

where: 

• x(t) is the particle position at the time i; 

• xq and vo are the initial position and velocity of 
the particle; 

• Ze and E are its charge and its (constant) energy; 

• B(x, t) is the magnetic field, which in principle 
can vary both with respect to the position and to 
the time. 

If it is possible to write down an explicit functional 
form for -B(x, t), an attempt to solve the equation [2] can 
be made. On the contrary, especially when the variation 
with the time cannot be easily summarized with an 
analytical formula, a numerical solution is unavoidably. 

The equation [2] explicitly shows the perturbation in- 
duced by the magnetic field on the straight trajectory 
(x(t) = xo + voi). This suggests an iterative method to 
determine the solution, which can be expressed as the 



where Ax 0(B ,+i)(t) = x 0(B ,+i)(t) - (x + v t) is 
the displacement from the unperturbed trajectory at the 
time t. Being content with the first-order approximation, 
it can be obtained: 

Zee 2 /"* [ T 

Ax(i) ~ — — v x I At I da B(x + v a, a) 
& Jo Jo 



Ax(i) ~ — v x T B (£;x ,v ) 



(3) 



where I&(t; xo, vo) is the integral of the magnetic field 
along the straight trajectory, whose value depends only 
on the time of the motion (t) and on its initial conditions 
(x and v ). 

Since the phenomenon studied concerns ultra- 
relativistic particles and fixing for a moment the initial 
position and the final time, the equation [5] becomes: 
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Ax ~ — v x 2b (v ) 
hi 

In short, on a first approximation the displacement de- 
pends only on the ratio charge-to-energy of the primary 
and on the initial direction of its ultrarelativistic motion 
(versor v )Q 

Let us consider only the lowest order of the geomag- 
netic field multipoles-expansion, i.e. the dipole term: 

3(b ■ x)x - x 2 h 

B(x) = — 

where b has intensity b ~ 8.1 ■ 10 27 T m 3 and the south 
magnetic pole is supposed to have coordinates 78.3° 
South, 111.0° East. By setting Vo||xo (vertical direction 
approximation) and integrating from YangBaJing to a 
distance ~ 60 Earth's radii, the equation Q] is immedi- 
ately obtained. 
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x(t) = X 0(B 0)(t) + X ( Bl )(t) + . . . 

where xo(B°)(t) = xo+vq£ is the unperturbed (straight) 
trajectory and for the higher orders holds: 

Ax ( B i+i)(t) = 



Zee 2 f , f T . dx ( B i) 



E 



Jo 



dr / da- 



do 



B(x 0(s ,),a) 



i = 0,l,. 

3 I.e. what has been done in the main part of this paper. 



4 The second dependence is not trivial at all. Because of the 
difference in the field-integrals, two identical particles having the same 
energy and starting together from the Moon can drift from the straight 
trajectory differently according to their different initial directions. 



